{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bikeshare数据集上的回归分析\n",
    "casual、registered和cnt三个特征均为要预测的y，作业里只需对cnt进行预测\n",
    "要求：最小二乘、岭回归、Lasso，比较结果说明原因\n",
    "\n",
    "思路：\n",
    "1.读入数据\n",
    "2.准备数据：0.8训练数据，0.2测试数据\n",
    "3.训练和评估"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 数据读取及基本处理\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "#模型\n",
    "from sklearn.linear_model import LinearRegression, RidgeCV, LassoCV, ElasticNetCV\n",
    "\n",
    "#模型评估\n",
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "#可视化\n",
    "import seaborn as sn\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# pandas display data frames as tables\n",
    "from IPython.display import display, HTML"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.读入数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>instant</th>\n",
       "      <th>season_1</th>\n",
       "      <th>season_2</th>\n",
       "      <th>season_3</th>\n",
       "      <th>season_4</th>\n",
       "      <th>mnth_1</th>\n",
       "      <th>mnth_2</th>\n",
       "      <th>mnth_3</th>\n",
       "      <th>mnth_4</th>\n",
       "      <th>mnth_5</th>\n",
       "      <th>...</th>\n",
       "      <th>weekday_5</th>\n",
       "      <th>weekday_6</th>\n",
       "      <th>temp</th>\n",
       "      <th>atemp</th>\n",
       "      <th>hum</th>\n",
       "      <th>windspeed</th>\n",
       "      <th>holiday</th>\n",
       "      <th>workingday</th>\n",
       "      <th>yr</th>\n",
       "      <th>cnt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0.355170</td>\n",
       "      <td>0.373517</td>\n",
       "      <td>0.828620</td>\n",
       "      <td>0.284606</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>985</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.379232</td>\n",
       "      <td>0.360541</td>\n",
       "      <td>0.715771</td>\n",
       "      <td>0.466215</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>801</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.171000</td>\n",
       "      <td>0.144830</td>\n",
       "      <td>0.449638</td>\n",
       "      <td>0.465740</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1349</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.175530</td>\n",
       "      <td>0.174649</td>\n",
       "      <td>0.607131</td>\n",
       "      <td>0.284297</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1562</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.209120</td>\n",
       "      <td>0.197158</td>\n",
       "      <td>0.449313</td>\n",
       "      <td>0.339143</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1600</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>6</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.180948</td>\n",
       "      <td>0.202329</td>\n",
       "      <td>0.532916</td>\n",
       "      <td>0.138482</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1606</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>7</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0.171197</td>\n",
       "      <td>0.170340</td>\n",
       "      <td>0.512798</td>\n",
       "      <td>0.301676</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1510</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>8</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0.131919</td>\n",
       "      <td>0.109191</td>\n",
       "      <td>0.550985</td>\n",
       "      <td>0.503869</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>959</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>9</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.098690</td>\n",
       "      <td>0.048706</td>\n",
       "      <td>0.446444</td>\n",
       "      <td>0.700017</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>822</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>10</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.114266</td>\n",
       "      <td>0.094271</td>\n",
       "      <td>0.496573</td>\n",
       "      <td>0.414115</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1321</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>10 rows × 35 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   instant  season_1  season_2  season_3  season_4  mnth_1  mnth_2  mnth_3  \\\n",
       "0        1         1         0         0         0       1       0       0   \n",
       "1        2         1         0         0         0       1       0       0   \n",
       "2        3         1         0         0         0       1       0       0   \n",
       "3        4         1         0         0         0       1       0       0   \n",
       "4        5         1         0         0         0       1       0       0   \n",
       "5        6         1         0         0         0       1       0       0   \n",
       "6        7         1         0         0         0       1       0       0   \n",
       "7        8         1         0         0         0       1       0       0   \n",
       "8        9         1         0         0         0       1       0       0   \n",
       "9       10         1         0         0         0       1       0       0   \n",
       "\n",
       "   mnth_4  mnth_5  ...  weekday_5  weekday_6      temp     atemp       hum  \\\n",
       "0       0       0  ...          0          1  0.355170  0.373517  0.828620   \n",
       "1       0       0  ...          0          0  0.379232  0.360541  0.715771   \n",
       "2       0       0  ...          0          0  0.171000  0.144830  0.449638   \n",
       "3       0       0  ...          0          0  0.175530  0.174649  0.607131   \n",
       "4       0       0  ...          0          0  0.209120  0.197158  0.449313   \n",
       "5       0       0  ...          0          0  0.180948  0.202329  0.532916   \n",
       "6       0       0  ...          1          0  0.171197  0.170340  0.512798   \n",
       "7       0       0  ...          0          1  0.131919  0.109191  0.550985   \n",
       "8       0       0  ...          0          0  0.098690  0.048706  0.446444   \n",
       "9       0       0  ...          0          0  0.114266  0.094271  0.496573   \n",
       "\n",
       "   windspeed  holiday  workingday  yr   cnt  \n",
       "0   0.284606        0           0   0   985  \n",
       "1   0.466215        0           0   0   801  \n",
       "2   0.465740        0           1   0  1349  \n",
       "3   0.284297        0           1   0  1562  \n",
       "4   0.339143        0           1   0  1600  \n",
       "5   0.138482        0           1   0  1606  \n",
       "6   0.301676        0           1   0  1510  \n",
       "7   0.503869        0           0   0   959  \n",
       "8   0.700017        0           0   0   822  \n",
       "9   0.414115        0           1   0  1321  \n",
       "\n",
       "[10 rows x 35 columns]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 读入数据，读取特征工程处理过后的数据（one-hot编码）\n",
    "data = pd.read_csv(\"FE_day.csv\")\n",
    "data.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 731 entries, 0 to 730\n",
      "Data columns (total 35 columns):\n",
      "instant         731 non-null int64\n",
      "season_1        731 non-null int64\n",
      "season_2        731 non-null int64\n",
      "season_3        731 non-null int64\n",
      "season_4        731 non-null int64\n",
      "mnth_1          731 non-null int64\n",
      "mnth_2          731 non-null int64\n",
      "mnth_3          731 non-null int64\n",
      "mnth_4          731 non-null int64\n",
      "mnth_5          731 non-null int64\n",
      "mnth_6          731 non-null int64\n",
      "mnth_7          731 non-null int64\n",
      "mnth_8          731 non-null int64\n",
      "mnth_9          731 non-null int64\n",
      "mnth_10         731 non-null int64\n",
      "mnth_11         731 non-null int64\n",
      "mnth_12         731 non-null int64\n",
      "weathersit_1    731 non-null int64\n",
      "weathersit_2    731 non-null int64\n",
      "weathersit_3    731 non-null int64\n",
      "weekday_0       731 non-null int64\n",
      "weekday_1       731 non-null int64\n",
      "weekday_2       731 non-null int64\n",
      "weekday_3       731 non-null int64\n",
      "weekday_4       731 non-null int64\n",
      "weekday_5       731 non-null int64\n",
      "weekday_6       731 non-null int64\n",
      "temp            731 non-null float64\n",
      "atemp           731 non-null float64\n",
      "hum             731 non-null float64\n",
      "windspeed       731 non-null float64\n",
      "holiday         731 non-null int64\n",
      "workingday      731 non-null int64\n",
      "yr              731 non-null int64\n",
      "cnt             731 non-null int64\n",
      "dtypes: float64(4), int64(31)\n",
      "memory usage: 200.0 KB\n"
     ]
    }
   ],
   "source": [
    "data.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "查看数据，没有数据缺失的情况"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.准备数据\n",
    "经过数据探索和特征工程(2. FE_BikeSharing)，没有数据缺失的情况，不需要处理。\n",
    "1.分离X,y\n",
    "2.分割训练数据、测试数据\n",
    "3.去掉ID‘instant’，ID不参与模型训练"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "#去掉X中的标签数据\n",
    "#drop函数默认删除行，删除列加axis=1\n",
    "y = data['cnt']\n",
    "X = data.drop(['cnt'],axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(584, 34)\n",
      "(147, 34)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/yy/.local/lib/python3.6/site-packages/sklearn/model_selection/_split.py:2069: FutureWarning: From version 0.21, test_size will always complement train_size unless both are specified.\n",
      "  FutureWarning)\n"
     ]
    }
   ],
   "source": [
    "# 用train_test_split 分割训练数据和测试数据\n",
    "from sklearn.model_selection import train_test_split\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, train_size = 0.8,random_state = 0)\n",
    "\n",
    "print(X_train.shape)\n",
    "print(X_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/yy/.local/lib/python3.6/site-packages/pandas/core/frame.py:3940: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  errors=errors)\n"
     ]
    }
   ],
   "source": [
    "#保存测试ID，用于结果提交\n",
    "testID = X_test['instant']\n",
    "\n",
    "#ID不参与预测\n",
    "#drop函数inplace = True，原数组名对应的内存值改变\n",
    "X_train.drop(['instant'], axis=1, inplace = True)\n",
    "X_test.drop(['instant'], axis=1, inplace = True)\n",
    "\n",
    "#保存特征名字以备后用（可视化）\n",
    "feat_names = X_train.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.模型训练和模型评估"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### （1） 最小二乘线性回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE on Training set : 767.8802034492847\n",
      "RMSE on Test set : 817.6251864872016\n"
     ]
    }
   ],
   "source": [
    "# 最小二乘没有超参数需要调优，直接用全体训练数据训练模型\n",
    "\n",
    "# 1. 生成学习器实例\n",
    "lr = LinearRegression()\n",
    "\n",
    "#2. 在训练集上训练学习器\n",
    "lr.fit(X_train, y_train)\n",
    "\n",
    "#3.用训练好的学习器对训练集/测试集进行预测\n",
    "y_train_pred = lr.predict(X_train)\n",
    "y_test_pred = lr.predict(X_test)\n",
    "\n",
    "rmse_train = np.sqrt(mean_squared_error(y_train,y_train_pred))\n",
    "rmse_test = np.sqrt(mean_squared_error(y_test,y_test_pred))\n",
    "print(\"RMSE on Training set :\", rmse_train)\n",
    "print(\"RMSE on Test set :\", rmse_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### （2）岭回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best alpha : 1.0\n",
      "cv of rmse : [ 806.14274372  806.13059312  806.011119    804.98946948  801.06246694\n",
      "  797.95313648  822.15888399 1078.04194503 1680.07803178 1878.71800678\n",
      " 1903.12118578 1905.61749753]\n",
      "RMSE on Training set : 754.0366623762044\n",
      "RMSE on Test set : 776.975360713381\n"
     ]
    }
   ],
   "source": [
    "# 1. 设置超参数搜索范围，生成学习器实例\n",
    "#RidgeCV多个alphas，得出多个对应最佳的超参数w,然后得到最佳的w及对应的alphas\n",
    "alphas = 10**np.linspace(-5, 6, 12)\n",
    "ridge = RidgeCV(alphas = alphas,store_cv_values=True )\n",
    "\n",
    "# 2. 用训练数据度模型进行训练\n",
    "ridge.fit(X_train, y_train)\n",
    "\n",
    "#通过交叉验证得到的最佳超参数alpha\n",
    "alpha = ridge.alpha_\n",
    "print(\"Best alpha :\", alpha)\n",
    "\n",
    "# 交叉验证估计的测试误差\n",
    "mse_cv = np.mean(ridge.cv_values_, axis = 0)\n",
    "rmse_cv = np.sqrt(mse_cv)\n",
    "print(\"cv of rmse :\",rmse_cv)\n",
    "\n",
    "#训练上测试，训练误差\n",
    "y_train_pred = ridge.predict(X_train)\n",
    "rmse_train = np.sqrt(mean_squared_error(y_train,y_train_pred))\n",
    "\n",
    "y_test_pred = ridge.predict(X_test)\n",
    "rmse_test = np.sqrt(mean_squared_error(y_test,y_test_pred))\n",
    "\n",
    "print(\"RMSE on Training set :\", rmse_train)\n",
    "print(\"RMSE on Test set :\" ,rmse_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAEKCAYAAAC7c+rvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XucVXW9//HXe2YYQG4DMiICCipe0Ap1RDxamaaix6LOz2PWMcksu+g59judX6mdDmV1Tj36nexux/KCHss8Vj/JvERmdmOQQRQFNYdBBUTAGa5ynZnP74/9HdzQDAzInjV77/fz8diPvfZnfdf6fpaiH9b6ftdaigjMzMyyUJF1AmZmVr5chMzMLDMuQmZmlhkXITMzy4yLkJmZZcZFyMzMMuMiZGZmmXERMjOzzLgImZlZZqqyTqC3Gz58eIwdOzbrNMzMisq8efNejYjaPbVzEdqDsWPH0tDQkHUaZmZFRdKL3Wnny3FmZpYZFyEzM8uMi5CZmWXGRcjMzDLjImRmZplxETIzs8y4CJmZWWZ8n5CZWRFrbw+2t7ezvS1obWtnW1tueXtrO9s7lttyyx3rWnf8zmvXnr9Nbl3dYUN521F7vN/0DXERMjPrpVau38LjL65h/tK1zH9pDUtbNu8oJq2puLS2R8H6/8QZR7gImZmVgy3b21j48nrmv7SG+S/lis7L67YAUF1ZwfGjBnP6+OH0raqgT2UF1VUVVFVox3Kfytxy7pO/XEF1laiqeH05f11+2+rKCvrsaCskFfy4XYTMzHpYRLB87WYeT8Vm/ktrWfTyera1tQMwqqY/J40dxkfG1HDCoTVMOGQwfasqM866MFyEzMwKbPO2NhYsW8v8pWt3XF5bvWErAP36VPDm0TVcdvpYTjx0KCeMqeGgwf0yzrjnuAiZme1HEcGLzZuYv3QNj7+4lvlL1/DMig20pbGbsQcewOlHDufEQ2s44dChHH3wIPpUlu9E5YIVIUn9gN8DfVM/90TEdEm3AW8H1qWmH4qIJ5S7+Pgt4HxgU4o/nvY1DfjX1P7LETEjxU8CbgP6A/cDV0dESBoG/BQYC7wAXBQRa3bXh5nZvti4tZUnl75+WW3+0rW0vLYNgAHVlUw8tIZPvP0ITkhFZ9iA6owz7l0KeSa0FTgzIjZK6gP8UdIDad3/iYh7dml/HjA+fU4BbgROSQVlOlAHBDBP0syIWJPafBSYQ64ITQEeAK4BHo6Ir0q6Jv3+bFd9FOTozaykLV69kavvms/Cl9cTaYLakQcN5KxjDuLEw4ZywqE1jD9oEJUVhR/cL2YFK0IREcDG9LNP+uxuLuFU4Pa0Xb2kGkkjgTOAWRHRAiBpFjBF0u+AwRFRn+K3A+8hV4Smpu0AZgC/I1eEOu0jIlbsl4M2s7IQEVz386dY2rKZq88azwmHDmXimBqG9O+TdWpFp6AXIiVVSnoCWEWukMxJq74iaYGkGyT1TbFRwNK8zZel2O7iyzqJA4zIKyyvACP20IeZWbf9csEK5ixp4TNTjuZT7zyKtx9V6wK0jwpahCKiLSImAqOBSZKOB64FjgFOBoaRO0MpZA7B7s/A/oqkKyQ1SGpYvXp1gTIzs2L02tZWvvKrRRw/ajAXn3xo1ukUvR6ZkhERa4FHgCkRsSJytgK3ApNSs+XAmLzNRqfY7uKjO4kDrEyX8kjfq/bQx6753hQRdRFRV1tb2LuFzay4fPeRRlau38oX3328x3v2g4IVIUm1kmrScn/gbODZvOIgcmM4T6dNZgKXKmcysC5dUnsIOEfSUElDgXOAh9K69ZImp31dCtybt69paXnaLvHO+jAz26Om1Rv50R+auPCk0Zx02NCs0ykJhZwdNxKYIamSXLG7OyLuk/RbSbWAgCeAj6f295ObOt1Ibvr0ZQAR0SLpS8Dc1O76jkkKwCd5fYr2A+kD8FXgbkmXAy8CF+2uDzOzPYkIvvDLRfSrquSzU47JOp2SoYjCPfyuFNTV1UVDQ0PWaZhZxn698BWuuGMen79gApefPi7rdHo9SfMiom5P7cr3Nl0zs27asr2N6+9bxFEjBnLpqYdlnU5J8WN7zMz24L8ebWLZms38+KOnlPUjdgrB/zTNzHZjacsmvv+7Ri5480j+5ojhWadTclyEzMx248u/WkSFxOf+9tisUylJLkJmZl149C+reWjhSq4680hGDumfdTolyUXIzKwT21rb+eLMhYw98AA+8lbPhisUT0wwM+vErX9aQtOrr3HrZSeX7FtNewOfCZmZ7eKVdVv49sPP885jR/COow/KOp2S5iJkZraL/3jgGba3B/92wYSsUyl5LkJmZnnmNDVz7xMv8/G3Hc6hBx6QdTolz0XIzCxpbWtn+syFjKrpzyfOODLrdMqCi5CZWXLnnJd49pUNfP6CY+lf7ckIPcFFyMwMeHXjVv7z18/x1vHDOfe4g7NOp2y4CJmZAV9/8Dk2bWtj+ruOI/eKMusJLkJmVvaeWLqWnzYs5cOnj+PIgwZmnU5ZcREys7LW3h5Mv/dpDhrUl38805MRepqLkJmVtf+Zt5Qnl63juvOPZVC/PlmnU3ZchMysbK3btJ2vPfgcJ48dytSJh2SdTllyETKzsvWNWc+xdtM2vvju4z0ZISMuQmZWlha9vJ476l/kksmHMeGQwVmnU7ZchMys7EQEX5i5kJoDqvnns4/KOp2y5iJkZmVn5pMv89gLLXzm3KOpOaA663TKmouQmZWVjVtb+cqvnuEto4dwUd2YrNMpe36pnZmVle88/DyrNmzlpkvrqKjwZISsFexMSFI/SY9JelLSQklfTPFxkuZIapT0U0nVKd43/W5M68fm7evaFH9O0rl58Skp1ijpmrz4XvdhZqWvcdVGbv7jEi6qG83EMTVZp2MU9nLcVuDMiHgLMBGYImky8DXghog4ElgDXJ7aXw6sSfEbUjskTQAuBo4DpgDfl1QpqRL4HnAeMAF4f2rL3vZhZqUvIvjiLxfSv7qSz0w5Jut0LClYEYqcjelnn/QJ4EzgnhSfAbwnLU9Nv0nrz1Ju4v5U4K6I2BoRS4BGYFL6NEZEU0RsA+4CpqZt9rYPMytxDy1cyR+ef5VPn30Uwwf2zTodSwo6MSGdsTwBrAJmAYuBtRHRmposA0al5VHAUoC0fh1wYH58l226ih+4D32YWQnbvK2NL923iGMOHsQlkw/LOh3LU9AiFBFtETERGE3uzKUozoElXSGpQVLD6tWrs07HzN6gGx9dzPK1m/niu4+jqtKTgnuTHvm3ERFrgUeAU4EaSR2z8kYDy9PycmAMQFo/BGjOj++yTVfx5n3oY9d8b4qIuoioq62t3cejNrPe4KXmTfzg0cW8+y2HcMrhvvDR2xRydlytpJq03B84G3iGXDG6MDWbBtyblmem36T1v42ISPGL08y2ccB44DFgLjA+zYSrJjd5YWbaZm/7MLMSdf19i6iqENedf2zWqVgnCnmf0EhgRprFVgHcHRH3SVoE3CXpy8B84ObU/mbgDkmNQAu5okJELJR0N7AIaAWujIg2AElXAQ8BlcAtEbEw7euze9OHmZWmR55bxW+eWck15x3DwUP6ZZ2OdUI+Edi9urq6aGhoyDoNM9tLW1vbmPLNPyDgwU+9jeoqjwX1JEnzIqJuT+38xAQzK0k3/3EJS159jds/PMkFqBfzvxkzKzkr1m3mOw83cu5xI3jbUZ5c1Ju5CJlZyfnKr56hPYJ//dsJe25smXIRMrOSMntxM/ctWMEnzjiCMcMOyDod2wMXITMrGdvb2vnCzIWMHtqfj7/9iKzTsW5wETKzknHH7Bd5buUG/u2CCfTrU5l1OtYNLkJmVhK2t7Xzzd/8hbcdVcvZE0ZknY51k4uQmZWEBcvWsX5LK+8/eQx+OH7xcBEys5JQ35R7DKSfD1dcXITMrCTUNzVzzMGDGDagOutUbC+4CJlZ0dvW2k7DC2uY7LOgouMiZGZF78lla9m8vY1Tj3ARKjYuQmZW9OoXNyPBKeOGZZ2K7SUXITMrerObmjn24MHUHODxoGLjImRmRW1raxvzXvR4ULFyETKzovbES2vZ2tru8aAi5SJkZkVtdlNuPGiSx4OKkouQmRW1+qZmjjtkMEP698k6FdsHLkJmVrS2bG/j8ZfWcqrHg4qWi5CZFa3HX1rDNo8HFTUXITMrWvVNLVQI6sZ6PKhYuQiZWdGqX9zMm0YNYXA/jwcVKxchMytKm7e1MX+p7w8qdi5CZlaUHn9pDdvbgskeDypqBStCksZIekTSIkkLJV2d4l+QtFzSE+lzft4210pqlPScpHPz4lNSrFHSNXnxcZLmpPhPJVWneN/0uzGtH7unPsysuMxe3ExlhTjZ40FFrZBnQq3ApyNiAjAZuFLShLTuhoiYmD73A6R1FwPHAVOA70uqlFQJfA84D5gAvD9vP19L+zoSWANcnuKXA2tS/IbUrss+CvePwMwKpb4pNx40sG9V1qnYG1CwIhQRKyLi8bS8AXgGGLWbTaYCd0XE1ohYAjQCk9KnMSKaImIbcBcwVbn3954J3JO2nwG8J29fM9LyPcBZqX1XfZhZEdm0rZUnl6311OwS0CNjQuly2AnAnBS6StICSbdIGppio4CleZstS7Gu4gcCayOidZf4TvtK69el9l3ty8yKSMMLaTzIkxKKXsGLkKSBwM+AT0XEeuBG4AhgIrAC+M9C57C3JF0hqUFSw+rVq7NOx8x2Ud/UTFWFqDts6J4bW69W0CIkqQ+5AnRnRPwcICJWRkRbRLQDP+T1y2HLgTF5m49Osa7izUCNpKpd4jvtK60fktp3ta+dRMRNEVEXEXW1tbX7cuhmVkCzm5p5y5gaBng8qOgVcnacgJuBZyLiG3nxkXnN3gs8nZZnAhenmW3jgPHAY8BcYHyaCVdNbmLBzIgI4BHgwrT9NODevH1NS8sXAr9N7bvqw8yKxMatrSxYto7Jh3tWXCko5F8jTgM+CDwl6YkUu47c7LaJQAAvAB8DiIiFku4GFpGbWXdlRLQBSLoKeAioBG6JiIVpf58F7pL0ZWA+uaJH+r5DUiPQQq5w7bYPMysODS+00NYenHr48KxTsf1AuRME60pdXV00NDRknYaZJf/xwDPc8sclLJh+Lv2rfYdFbyVpXkTU7amdn5hgZkWlvqmFiWNqXIBKhIuQmRWNDVu28/TydX5/UAlxETKzojE3jQf5/qDS4SJkZkWjvqmF6soKTvT9QSXDRcjMisbsxc2ccGgN/fp4PKhUuAiZWVFYt3k7C19e50txJcZFyMyKwtwlLbQHfmhpiXERMrOiMLupmb5VFUwcU5N1KrYfuQiZWVGob2rmxEOHejyoxLgImVmvt3bTNhatWO9LcSXIRcjMer05S1qIwJMSSlC3i5Ck0yVdlpZr01OozcwKrr6pmX59KnjLmCFZp2L7WbeKkKTp5J5YfW0K9QH+u1BJmZnlm724mbrDhtG3yuNBpaa7Z0LvBd4NvAYQES8DgwqVlJlZhzWvbePZVzb4/UElqrtFaFt6KVwASBpQuJTMzF43Z0kz4PuDSlV3i9Ddkv6L3Ou0Pwr8htyruc3MCmr24mb696nkzaN9f1Ap6tabVSPi/0o6G1gPHA38W0TMKmhmZmbkHlpaN3YofSo9mbcUdasIpctvv42IWZKOBo6W1Ccithc2PTMrZ80bt/Lcyg1MPeGQrFOxAunuXy1+D/SVNAp4EPggcFuhkjIzg9xZEPj+oFLW3SKkiNgE/B1wY0T8PXBc4dIyM8vdHzSgupI3jfL9QaWq20VI0qnAPwC/SjFP2Dezgprd1MzJ44Z5PKiEdfff7NXANcDPI2JhelrCbwuXlpmVu9UbttK4aqMvxZW4bk1MADYB7cD7JV0CiHTPkJlZIdQ3pfuDXIRKWneL0J3AvwBPkytGZmYFNbupmYF9qzjukMFZp2IF1N3Lcasj4pcRsSQiXuz47G4DSWMkPSJpkaSFkq5O8WGSZkl6Pn0PTXFJ+rakRkkLJJ2Yt69pqf3zkqblxU+S9FTa5tuStK99mFnvUt/UzKRxw6jyeFBJ6+6/3emSfiTp/ZL+ruOzh21agU9HxARgMnClpAnkxpYejojxwMPpN8B5wPj0uQK4EXIFBZgOnAJMSrkMTdvcCHw0b7spKb5XfZhZ77Jy/RaaVr/mS3FloLtF6DJgIrn/yb8rfS7Y3QYRsSIiHk/LG4BngFHAVGBGajYDeE9angrcHjn15B4RNBI4F5gVES0RsQaYBUxJ6wZHRH16rt3tu+xrb/ows16kYzzIkxJKX3fHhE6OiKP3tRNJY4ETgDnAiIhYkVa9AoxIy6OApXmbLUux3cWXdRJnH/pYgZn1GvVNzQzuV8UEjweVvO6eCf05XUrba5IGAj8DPhUR6/PX5T+Zu1D2pQ9JV0hqkNSwevXqAmVmZl2ZvbiZSeMOpLJCWadiBdbdIjQZeELSc2lA/ylJC/a0kaQ+5ArQnRHx8xRe2XEJLH2vSvHlwJi8zUen2O7iozuJ70sfO4mImyKiLiLqamtr93SYZrYfrVi3mReaN/n9QWWiu0VoCrnB/HN4fTzoXbvbIM1Uuxl4JiK+kbdqJtAxw20acG9e/NI0g20ysC5dUnsIOEfS0DQh4RzgobRuvaTJqa9Ld9nX3vRhZr3EjvuD/P6gstDdVznsdjp2F04j96DTpyQ9kWLXAV8l936iy4EXgYvSuvuB84FGcjfHXpb6bpH0JWBuand9RLSk5U+Se5Bqf+CB9GFv+zCz3mP24maG9O/DsQd7PKgcdHdiwl6LiD+Se7JCZ87qpH0AV3axr1uAWzqJNwDHdxJv3ts+zKx3qG9q4ZRxw6jweFBZ8F1gZtZrLF+7mZdaNvlSXBlxETKzXqN+se8PKjcuQmbWa8xuamboAX04esSgrFOxHuIiZGa9xuzFzUw+/ECPB5URFyEz6xWWtmxi+drNvhRXZlyEzKxXmO37g8qSi5CZ9Qr1i5s5cEA14w8amHUq1oNchMwscxFBfVNuPCi9FszKhIuQmWXupZZNvLxuC5N9Ka7suAiZWeZ2PC/ODy0tOy5CZpa52YubqR3UlyNqPR5UblyEzCxTEcFsjweVLRchM8vUC82bWLl+q98fVKZchMwsU7MXd4wHeVJCOXIRMrNMzW5q5qBBfRk3fEDWqVgGXITMLDMd9wedeoTHg8qVi5CZZWbx6tdYvWGrL8WVMRchM8tMx/1Bfmhp+XIRMrPMzG5qZuSQfhx24AFZp2IZcREys0xEBHN8f1DZcxEys0w0rtrIqxu3eTyozLkImVkm/P4gAxchM8vI7MXNjKrpz+ih/bNOxTLkImRmPa69PZizpMXjQVa4IiTpFkmrJD2dF/uCpOWSnkif8/PWXSupUdJzks7Ni09JsUZJ1+TFx0mak+I/lVSd4n3T78a0fuye+jCznvWXVRtoeW2bL8VZQc+EbgOmdBK/ISImps/9AJImABcDx6Vtvi+pUlIl8D3gPGAC8P7UFuBraV9HAmuAy1P8cmBNit+Q2nXZx34+ZjPrhvrFHfcH+aGl5a5gRSgifg+0dLP5VOCuiNgaEUuARmBS+jRGRFNEbAPuAqYqd/5+JnBP2n4G8J68fc1Iy/cAZ6X2XfVhZj1sdlMzY4b1Z/RQ3x9U7rIYE7pK0oJ0uW5oio0Clua1WZZiXcUPBNZGROsu8Z32ldavS+272peZ9aAd40HjfCnOer4I3QgcAUwEVgD/2cP9d4ukKyQ1SGpYvXp11umYlZRnX9nA2k3bPR5kQA8XoYhYGRFtEdEO/JDXL4ctB8bkNR2dYl3Fm4EaSVW7xHfaV1o/JLXval+d5XlTRNRFRF1tbe2+HKqZdWG2nxdneXq0CEkamffzvUDHzLmZwMVpZts4YDzwGDAXGJ9mwlWTm1gwMyICeAS4MG0/Dbg3b1/T0vKFwG9T+676MLMeNHtxM4cdeACH1Pj+IIOqPTfZN5J+ApwBDJe0DJgOnCFpIhDAC8DHACJioaS7gUVAK3BlRLSl/VwFPARUArdExMLUxWeBuyR9GZgP3JziNwN3SGokNzHi4j31YWY9o609eGxJM+e/aeSeG1tZUO4kwbpSV1cXDQ0NWadhVhKeXr6OC77zR7518USmTvS8oFImaV5E1O2pnZ+YYGY9xu8Psl25CJlZj5m9uJnDhw9gxOB+WadivYSLkJn1iNa2dh5b0sIpPguyPC5CZtYjFq1Yz4atrb4/yHbiImRmPWK2nxdnnXARMrMeUd/UzBG1AzhokMeD7HUuQmZWcK1t7cx9YY0vxdlfcREys4J7avk6Nm5t5dTDh2edivUyLkJmVnD1Tbm3upzi8SDbhYuQmRXc7KZmjhoxkOED+2adivUyLkJmVlDb29ppeKHFT0mwTrkImVlBLVi2jk3b2jjVRcg64SJkZgXV8bw4PynBOuMiZGYFVd/UzDEHD2LYgOqsU7FeyEXIzApmW2s7DS+s8XiQdclFyMwK5slla9m8vc03qVqXXITMrGD+3NiMBKeM8/1B1jkXITMriGdfWc8P/9DEyWOHUXOAx4Oscy5CZrbfrVy/hQ/fOpcBfSv51sUTs07HerGqrBMws9KyaVsrl8+Yy9rN27n7Y6cyckj/rFOyXsxnQma237S1B//0k/ksenk93/vAiRw/akjWKVkv5zMhM9tvvnTfIn7zzCq+NPU43nHMQVmnY0XAZ0Jmtl/c+qcl3PbnF7j89HF88NSxWadjRcJFyMzesFmLVnL9fYs497gRXHf+sVmnY0WkYEVI0i2SVkl6Oi82TNIsSc+n76EpLknfltQoaYGkE/O2mZbaPy9pWl78JElPpW2+LUn72oeZ7bunlq3jn34ynzePGsI333cClRXKOiUrIoU8E7oNmLJL7Brg4YgYDzycfgOcB4xPnyuAGyFXUIDpwCnAJGB6R1FJbT6at92UfenDzPbd8rWb+fCMuQwbUM0Pp9XRv7oy65SsyBSsCEXE74GWXcJTgRlpeQbwnrz47ZFTD9RIGgmcC8yKiJaIWAPMAqakdYMjoj4iArh9l33tTR9mtg/Wb9nOh2+dy5btbdx62ckcNKhf1ilZEerpMaEREbEiLb8CjEjLo4Clee2Wpdju4ss6ie9LH2a2l7a3tXPlnY+zePVGfnDJSRw1YlDWKVmRymxiQjqDid7Yh6QrJDVIali9enUBMjMrXhHBv937NH94/lX+/b1v4rQjh2edkhWxni5CKzsugaXvVSm+HBiT1250iu0uPrqT+L708Vci4qaIqIuIutra2r06QLNS94NHm/jJY0u58h1HcNHJY/a8gdlu9HQRmgl0zHCbBtybF780zWCbDKxLl9QeAs6RNDRNSDgHeCitWy9pcpoVd+ku+9qbPsysm+5b8DJfe/BZ3vWWQ/j02UdnnY6VgII9MUHST4AzgOGSlpGb5fZV4G5JlwMvAhel5vcD5wONwCbgMoCIaJH0JWBuand9RHRMdvgkuRl4/YEH0oe97cPMumfei2v457ufpO6woXz9wjdT4anYth8oN2xiXamrq4uGhoas0zDL1IvNr/He7/+Zwf2q+PknT/Orum2PJM2LiLo9tfMTE8xst9Zu2sZlt82lPYJbL5vkAmT7lYuQmXVpa2sbV9wxj2Utm/nhpXWMGz4g65SsxPgp2mbWqYjg2p89xWNLWvjWxRM5eaxf0W37n8+EzKxT3/zN8/x8/nI+ffZRTJ3o+7qtMFyEzOyv/GzeMr718PNceNJorjrzyKzTsRLmImRmO5m9uJlrfr6AvzniQP79vW8iPaDerCBchMxsh8ZVG/nYHQ0cduAAbrzkJKqr/L8IKyz/CTMzAF7duJXLbnuM6qoKbv3QyQzp3yfrlKwMeHacmbFlexsfvb2BVeu38tOPncqYYQdknZKVCRchszLX3h78891P8MTStdz4DycycUxN1ilZGfHlOLMy97WHnuX+p17huvOOZcrxfs+j9SwXIbMy9uM5L/FfjzZxyeRD+chbx2WdjpUhFyGzMvW751bx+Xuf5oyja/nCu47zVGzLhIuQWRl6ZsV6rvrxfI4aMYjvfuBEqir9vwLLhv/kmZWZleu38OHb5jKgbyW3fKiOgX09P8my4z99ZiUsItiwtZXVG7bu+Pzg0cWs27yd//n4qYwc0j/rFK3MuQiZFaEt29tyRWXj1p0KTP7vV9Py1tb2nbatrqrgB5ecyHGHDMkoe7PXuQgVyNwXWvjB7xZnnUanujf+vOdGUq5V7lu577TMjnXKa/P6b/K36WQf7LKuQqKyQlRKVFaKqgpRWVGRvrXj+/XlvHWVedvu+L3zthV5+6iqqKCy4vU+d/SdvyxRUUEnsX0f3G9ta6fltW2syismr3ZRZDZsae10HwcOqKZ2UF9qB/Xl8OEDdizXDurL8IG575FD+jGon5+GYL2Di1CBbNnexsoNW7JO4690523u3WpD7lJPR/sg0ncuHqlR/u+d2kVHX52s25HD67/bI2hr3/nT2t47X02/U5FKhSm/SHUUwx3rJdZv2U7za9s6/Wc/sG9VrpAM7MuxBw/mbeP77vidX2SGDaimjycYWJFxESqQt46v5a3ja7NOo6RFBO0Bre3tO4pSW1vQFrHT753W5xWw3HfX27Z3tI+0HK/HWtsjFcadC+SO5R1t6SSWv8/cEwsG96/6q6JSO7AfwwdVc0C1/zO10uU/3Va0JFEpqKyozDoVM9tHPnc3M7PMuAiZmVlmXITMzCwzmRQhSS9IekrSE5IaUmyYpFmSnk/fQ1Nckr4tqVHSAkkn5u1nWmr/vKRpefGT0v4b07baXR9mZpaNLM+E3hEREyOiLv2+Bng4IsYDD6ffAOcB49PnCuBGyBUUYDpwCjAJmJ5XVG4EPpq33ZQ99GFmZhnoTZfjpgIz0vIM4D158dsjpx6okTQSOBeYFREtEbEGmAVMSesGR0R95G5kuX2XfXXWh5mZZSCrIhTAryXNk3RFio2IiBVp+RVgRFoeBSzN23ZZiu0uvqyT+O76MDOzDGR1n9DpEbFc0kHALEnP5q+MiJBU0Nvhd9dHKoxXABx66KGFTMPMrKxlUoQiYnn6XiXpF+TGdFZKGhkRK9IltVWp+XJgTN7YuoofAAAGSUlEQVTmo1NsOXDGLvHfpfjoTtqzmz52ze8m4CYASaslvbivx9rDhgOvZp1EgZTysUFpH5+PrXi9keM7rDuNerwISRoAVETEhrR8DnA9MBOYBnw1fd+bNpkJXCXpLnKTENalIvIQ8O95kxHOAa6NiBZJ6yVNBuYAlwLfydtXZ310KSKK5tk7khryJnqUlFI+Nijt4/OxFa+eOL4szoRGAL9Is6argB9HxIOS5gJ3S7oceBG4KLW/HzgfaAQ2AZcBpGLzJWBuand9RLSk5U8CtwH9gQfSB3LFp7M+zMwsA4ruPDLZikIp/62slI8NSvv4fGzFqyeOrzdN0bY37qasEyigUj42KO3j87EVr4Ifn8+EzMwsMz4TMjOzzLgIlShJn5YUkoZnncv+Iunrkp5NzxD8haSarHN6oyRNkfRces5hST1GStIYSY9IWiRpoaSrs85pf5NUKWm+pPuyzmV/klQj6Z7039szkk4tVF8uQiVI0hhyU9ZfyjqX/WwWcHxEvBn4C3Btxvm8IZIqge+Rez7iBOD9kiZkm9V+1Qp8OiImAJOBK0vs+ACuBp7JOokC+BbwYEQcA7yFAh6ji1BpugH4DLnHI5WMiPh1RLSmn/XsfFNyMZoENEZEU0RsA+4i93zDkhARKyLi8bS8gdz/yEbtfqviIWk08LfAj7LOZX+SNAR4G3AzQERsi4i1herPRajESJoKLI+IJ7POpcA+zOv3fxWrrp5/WHIkjQVOIHcDean4Jrm/7LVnnch+Ng5YDdyaLjX+KD1YoCCyenacvQGSfgMc3MmqzwHXkbsUV5R2d2wRcW9q8zlyl3ru7MncbN9IGgj8DPhURKzPOp/9QdIFwKqImCfpjKzz2c+qgBOBf4yIOZK+Re61N58vVGdWZCLinZ3FJb2J3N9inkxPpBgNPC5pUkS80oMp7rOujq2DpA8BFwBnRfHfX9DVcxFLhqQ+5ArQnRHx86zz2Y9OA94t6XygHzBY0n9HxCUZ57U/LAOWRUTHWes9FPDda75PqIRJegGoi4iSeMCipCnAN4C3R8TqrPN5oyRVkZtgcRa54jMX+EBELMw0sf0kvdF4BtASEZ/KOp9CSWdC/xIRF2Sdy/4i6Q/ARyLiOUlfAAZExP8pRF8+E7Ji8l2gL7nXfwDUR8THs01p30VEq6SrgIeASuCWUilAyWnAB4GnJD2RYtdFxP0Z5mTd84/AnZKqgSbSMzsLwWdCZmaWGc+OMzOzzLgImZlZZlyEzMwsMy5CZmaWGRchMzPLjIuQWYFI2vgGt79H0uF7aPM7Sbt982V32uzSvlbSg91tb/ZGuAiZ9UKSjgMqI6Kpp/tONwKvkHRaT/dt5cdFyKzAlPN1SU9LekrS+1K8QtL30ztbZkm6X9KFabN/AO7N28eNkhrSe3m+2EU/GyXdkNo8LKk2b/XfS3pM0l8kvTW1HyvpD5IeT5+/yWv//1IOZgXlImRWeH8HTCT3XpZ3Al+XNDLFx5J7l9AHgfwXh50GzMv7/bmIqAPeDLxd0ps76WcA0BARxwGPAtPz1lVFxCTgU3nxVcDZEXEi8D7g23ntG4C37v2hmu0dP7bHrPBOB34SEW3ASkmPAien+P9ERDvwiqRH8rYZSe5x+h0uknQFuf9mR5IrXAt26acd+Gla/m8g/4GhHcvzyBU+gD7AdyVNBNqAo/LarwIO2cvjNNtrLkJmvdNmck9nRtI44F+AkyNijaTbOtbtQf4zubam7zZe/+/+fwMryZ2hVQBb8tr3SzmYFZQvx5kV3h+A90mqTOM0bwMeA/4E/K80NjQCOCNvm2eAI9PyYOA1YF1qd14X/VQAHWNKHwD+uIe8hgAr0pnYB8k9RLXDUcDT3Tg2szfEZ0JmhfcLcuM9T5I7O/lMRLwi6WfkXuOwiNwbVh8H1qVtfkWuKP0mIp6UNB94NrX7Uxf9vAZMkvSv5C6nvW8PeX0f+JmkS4EH0/Yd3pFyMCsoP0XbLEOSBkbERkkHkjs7Oi0VqP7AI+l3Wzf3tTEiBu6nvH4PTI2INftjf2Zd8ZmQWbbuk1QDVANf6ngDbkRsljQdGAW81JMJpUuG33ABsp7gMyEzM8uMJyaYmVlmXITMzCwzLkJmZpYZFyEzM8uMi5CZmWXGRcjMzDLz/wHwikapOKW9BAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#均值\n",
    "mse_mean = np.mean(ridge.cv_values_, axis = 0)\n",
    "plt.plot(np.log10(alphas), mse_mean.reshape(len(alphas),1)) \n",
    "\n",
    "plt.xlabel('log(alpha)')\n",
    "plt.ylabel('mse')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### （3）Lasso"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best alpha : 2.336453166053946\n",
      "cv of rmse : 828.5550460506832\n",
      "RMSE on Training set : 754.2656275844785\n",
      "RMSE on Test set : 786.623058576428\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/yy/.local/lib/python3.6/site-packages/sklearn/model_selection/_split.py:1943: FutureWarning: You should specify a value for 'cv' instead of relying on the default value. The default value will change from 3 to 5 in version 0.22.\n",
      "  warnings.warn(CV_WARNING, FutureWarning)\n"
     ]
    }
   ],
   "source": [
    "#1. 生成学习器实例，LassoCV默认参数可自动确定alpha的搜素范围\n",
    "lasso = LassoCV()\n",
    "\n",
    "#2.模型训练\n",
    "lasso.fit(X_train, y_train)\n",
    "alphas2 = lasso.alpha_\n",
    "print(\"Best alpha :\" , alphas2)\n",
    "\n",
    "#3. 模型性能：cv\n",
    "mse_cv = np.mean(lasso.mse_path_, axis = 1)\n",
    "rmse_cv = np.sqrt(mse_cv)\n",
    "print(\"cv of rmse :\", min(rmse_cv))\n",
    "\n",
    "#训练误差\n",
    "y_train_pred = lasso.predict(X_train)\n",
    "rmse_train = np.sqrt(mean_squared_error(y_train,y_train_pred))\n",
    "print(\"RMSE on Training set :\" ,rmse_train)\n",
    "\n",
    "#测试误差\n",
    "y_test_pred = lasso.predict(X_test)\n",
    "rmse_test = np.sqrt(mean_squared_error(y_test,y_test_pred))\n",
    "print(\"RMSE on Test set :\" ,rmse_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_test_pred = ridge.predict(X_test)\n",
    "\n",
    "#生成提交测试结果\n",
    "df = pd.DataFrame({\"instant\":testID, 'cnt':y_test_pred})\n",
    "df.to_csv('submission.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (1)OLS：系数的值非常大，OLS模型的性能并不好\n",
    "### (2)岭回归：相比OLS，岭回归模型增加了L2正则，系数值进行了收缩。 由于增加正则限制了模型复杂的，相比OLS模型，岭回归模型在测试集上的误差有所减小。\n",
    "### (3)Lasso：相比OLS，Lasso模型增加了L1正则，系数值进行了收缩，同时有些特征的系数为0。 在共享单车案例中，岭回归模型比Lasso模型性能稍微好一点。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
